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Abstract 

In this work we consider five different lattice models which exhibit continuous phase 
transitions into absorbing states. By measuring certain universal functions, which 
characterize the steady state as well as the dynamical scaling behavior, we present 
clear numerical evidence that all models belong to the universality class of directed 
percolation. Since the considered models are characterized by different interaction 
details the obtained universal scaling plots are an impressive manifestation of the 
universality of directed percolation. 
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1 Introduction 

In this work we consider the universality class of directed percolation (DP, 
see [1,2] for recent reviews). Because of its robustness and ubiquity, (includ- 
ing critical phenomena in physics, biology, epidemiology, as well as catalytic 
chemical reactions) directed percolation is recognized as the paradigm of non- 
equilibrium phase transitions into absorbing states. These so-called absorb- 
ing phase transitions arise from a competition of opposing processes, usu- 
ally creation and annihilation processes. The transition point separates an 
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active phase from an absorbing phase in which the dynamics is frozen. Anal- 
ogous to equilibrium critical phenomena, absorbing phase transitions can be 
grouped into different universality classes. All systems belonging to a given 
universality class share the same critical exponents, and certain scaling func- 
tions (e.g. equation of state, correlation functions, finite-size scaling functions, 
etc.) become identical near the critical point. According to the universality 
hypothesis of Janssen and Grassberger, short-range interacting models, ex- 
hibiting a continuous phase transition into a unique absorbing state, belong 
to the directed percolation universality class, if they are characterized by a 
one-component order parameter and no additional symmetries [3,4]. 

Similar to equilibrium critical phenomena, the universality of directed per- 
colation is understood by renormalization group treatments of an associated 
continuous field theory. The process of directed percolation might be repre- 
sented by the Langevin equation [3] 

d t n(x,t) = rn(x,t) — un 2 (x,t) + W 2 n(x, t) + r)(x,t) . (1) 



Here, n(x, t) corresponds to the density of active sites on a mesoscopic scale 
and r describes the distance to the critical point. Furthermore, 77 denotes the 
noise which accounts for fluctuations of n(x, t). According to the central limit 
theorem, r)(x,t) is a Gaussian random variable with zero mean and whose 
correlator is given by 

(r](x,t)r]{x!,t')) = Kn{x,t)S(x-x')S(t-t'). (2) 



Notice, that the latter equation ensures that the system is trapped in the ab- 
sorbing state n(x, t) = 0. Furthermore, higher order terms such as n(x, t) 3 , 
n(x, t) 4 or V 4 n(x, t) are irrelevant under renormalization group transforma- 
tions as long as u > 0. Negative values of u give rise to a first order phase 
transition whereas u = is associated with a tricritical point [5,6]. 

Above the upper critical dimension D c mean field theories apply and present 
instructive insight into the critical behavior [7]. Below D c , renormalization 
group techniques have to be applied to determine the critical exponents and 
the scaling functions (see [8,9] for recent reviews of the field theoretical treat- 
ment of directed percolation). In that case path integral formulations are more 
adequate than the Langevin equation. Stationary correlation functions as well 
as response functions can be determined by calculating path integrals with 
weight exp (— S), where the dynamic functional S describes the considered 
stochastic process. Up to higher irrelevant orders the dynamic functional as- 
sociated with directed percolation is given by [3,10,11,12] 



S[n,n] = J d D xdtn 



d t n — (r + V 2 )n — ( — n — u n ) n 



(3) 
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where n(x, t) denotes the response field conjugated to the Langevin noise 
field [13]. Remarkably, the above functional is well known from high energy 
physics and corresponds to the Lagrangian of Reggeon field theory [14]. Since 
DP represents the simplest realization of a nonequilibrium phase transition its 
field theory is often regarded as the nonequilibrium counterpart of the famous 
4 -theory of equilibrium [8]. Rescaling the fields 

n(x,t) = fj,s(x,t), n(x,t) = /i^ 1 s(x,t) (4) 



the functional S is invariant under the duality transformation (so-called ra- 
pidity reversal in Reggeon field theory) 

s(x,t) < — > —s(x,—t) (5) 



for /i 2 = 2u/k. Note that /i is a redundant variable from the renormalization 
group point of view [3,15]. The rapidity reversal (5) is the characteristic sym- 
metry of the universality class of directed percolation. It is worth to reempha- 
size that the rapidity reversal is obtained from a field theoretical treatment. 
Thus all models that belong to the universality class of directed percolation 
obey, at least asymptotically, the rapidity reversal after a corresponding coarse 
grained procedure [3]. In other words, the rapidity symmetry may not be repre- 
sented in the microscopic models. For example, bond directed percolation [1] 
and site directed percolation [16] obey the rapidity reversal microscopically 
whereas e.g. the contact process and the pair contact process do not. Still, the 
latter two models do so asymptotically and belong to the DP class. 

Due to the continuing improvement of computer hardware, high accurate nu- 
merical data of DP have become available in the last years, resulting in a 
fruitful and instructive interplay between numerical investigations and renor- 
malization group analyses. In particular, investigations of the scaling behavior 
of the equation of state and of the susceptibility [17,18,19], of finite-size scaling 
functions [20,21], and of dynamical scaling functions [18,22] yield an impres- 
sive agreement between field theoretical and numerical results. For example, 
the universal amplitude ratio of the susceptibility has been calculated via e- 
expansions [17]. In second order of e, the error of the field theory estimate 
is about 6% for D = 2 (see ref. [16] for a detailed discussion). It is instruc- 
tive to compare this result to the equilibrium situation. The corresponding 
<p -theory value [23] differs from the exact value of the two-dimensional Ising 
model [24,25] by roughly 115%. Thus, in contrast to the 4 -theory the DP field 
theory provides excellent numerical estimates of certain universal quantities. 

In our previous works we investigated the steady state critical scaling behav- 
ior of two one-dimensional models [26] and of one model in various dimen- 
sions [19], respectively. In this work we extend our investigations and consider 
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the steady state and dynamical scaling behavior of five different models in var- 
ious dimensions. All lattice models are expected to belong to the universality 
class of directed percolation. So far, most works focus on the determination 
of the critical exponents only, neglecting the determination of universal scal- 
ing functions. It turns out that checking the universality class it is often a 
more exact test to consider scaling functions rather than the values of the 
critical exponents. While for the latter ones the variations between different 
universality classes are often small, the scaling functions may differ signifi- 
cantly [27]. Thus the agreement of universal scaling functions provides not 
only additional but also independent and more convincing evidence in favor 
of the conjecture that the phase transitions of two models belong to the same 
universality class. Additionly to the critical exponents and scaling functions, 
universality classes are also characterized by certain amplitude combinations 
(see e.g. [28]). But these amplitude combinations are merely particular values 
of the scaling functions and will be neglected in this work. 



2 Lattice models of directed percolation 

In the following we consider various lattice models that belong to the univer- 
sality class of directed percolation. First, we revisit the contact process that is 
well known in the mathematical literature (see e.g. [29]). Second, we consider 
the Domany-Kinzel cellular automaton [30] which is very useful in order to 
perform numerical investigations of directed bond and directed site percola- 
tion. Third, we consider the pair contact process [31] that is characterized 
in contrast to the other models by infinitely many absorbing states. Unlike 
the first two models the universal scaling behavior of the pair contact pro- 
cess is still a matter of discussions in the literature. Furthermore, we briefly 
discuss the threshold transfer process [32] as well as the Ziff-Gulari-Barshad 
model [33]. The latter one mimics the catalysis of carbon monoxide oxidation. 

2. 1 Contact process 

The contact process (CP) is a continuous-time Markov process that is usu- 
ally defined on a D-dimensional simple cubic lattice (see [34] and references 
therein). A lattice site may be empty (n = 0) or occupied (n — 1) by a par- 
ticle and the dynamics is characterized by spontaneously occurring processes, 
taking place with certain transition rates. In numerical simulations the asyn- 
chronous update is realized by a random sequential update scheme: A particle 
on a randomly selected lattice site % is annihilated with rate one, whereas 
particle creation takes places on an empty neighboring site with rate XN/2D, 
i.e., 
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rii — 1 — j-> rij = , (6) 
= — > = 1 , (7) 

XN/2D ' W 

where A" denotes the number of occupied neighbors of ?v Note that the rates 
are defined as transition probabilities per time unit, i.e., they may be larger 
than one. Thus, rescaling the time will change the transition rates. In simula- 
tions a discrete time formulation of the contact process is performed. In that 
case a particle creation takes place at a randomly chosen neighbor site with 
probability p — A/(l + A) whereas particle annihilation occurs with proba- 
bility 1 — p = 1/(1 + A). In dynamical simulations, the time increment 1/A" a 
is associated with each attempted elementary update step, where iV a denotes 
the number of active sites. It is usual to present the critical value in terms of 
A c instead of p c . 

Similar to equilibrium phase transitions, it is often possible for absorbing 
phase transitions to apply an external field h that is conjugated to the order 
parameter, i.e., to the density of active sites p a . Being a conjugated field it 
has to destroy the absorbing phase, it has to be independent of the control 
parameter, and the corresponding linear response function has to diverge at 
the critical point 

X = M — oo. (8) 



In case of the CP, the conjugated field is implemented by a spontaneous cre- 
ation of particles, i.e., the external field creates a particle at an empty lattice 
site with rate h. Clearly spontaneous particle generation destroys the absorb- 
ing state and therefore the absorbing phase transition at all. Incorporating the 
conjugated field, a series of opportunities is offered to compare renormaliza- 
tion group results to those of numerical investigations. For example, simula- 
tions performed for non-zero field include the measurements of the equation of 
state [26], of the susceptibility [19], as well as of a modified finite-size scaling 
analysis appropriate for absorbing phase transitions [21,35]. 



2.2 Domany-Kinzel automaton 



An important 1 + 1-dimensional stochastic cellular automaton exhibiting di- 
rected percolation scaling behavior is the Domany-Kinzel (DK) automaton [30]. 
It is defined on a diagonal square lattice with a discrete time variable and 
evolves by parallel update according to the following rules: A site at time t is 
occupied with probability p 2 (pj if both (only one) backward sites (at time 
t — 1) are occupied. Otherwise the site remains empty. If both backward sites 
are empty a spontaneous particle creation takes place with probability p . 
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Similar to the contact process, the spontaneous particle creation destroys the 
absorbing phase (empty lattice) and corresponds therefore to the conjugated 
field h. 

It is straight forward to generalize the 1 + 1-dimensional Domany-Kinzel au- 
tomaton to higher dimensions (see e.g. [36,37,19]). In the following, we consider 
cellular automata onaD + 1-dimensional body centered cubic (bcc) lattice 
where the time corresponds to the [0, 0, ... , 0, 1] direction. A lattice site at 
time t is occupied with probability p if at least one of its 2 D backward neigh- 
boring sites (at time t — 1) is occupied. Otherwise the site remains empty. This 
parameter choice corresponds to the probabilities p x — p 2 — ■ ■ ■ — p 2 D = P, 
i.e., site-directed percolation (sDP) is considered. 



2.3 Pair contact process 

The pair contact process (PCP) was introduced by Jensen [31] and is one of 
the simplest models with infinitely many absorbing states showing a continu- 
ous phase transition. The process is defined on a D-dimensional cubic lattice 
and an asynchronous (random sequential) update scheme is applied. A lattice 
site may be either occupied (n — 1) or empty (n — 0). Pairs of adjacent occu- 
pied sites, linked by an active bond, annihilate each other with probability p 
otherwise an offspring is created at a neighboring site provided the target site 
is empty. The density of active bonds p a is the order parameter of a continuous 
phase transition from an active state to an inactive absorbing state without 
particle pairs. Similar to the contact process and to the Domany-Kinzel au- 
tomaton a spontaneous particle creation acts as a conjugated field [26]. Since 
isolated particles remain inactive, any configuration containing only isolated 
particles is absorbing. In case of the 1 + 1-dimensional pair contact process 
with L sites and periodic boundary conditions the number of absorbing states 
is asymptotically given by the golden mean N ~ (1/2 + y/E/2) L [38]. In the 
thermodynamic limit (L — > oo), the pair contact process is characterized by 
infinitely many absorbing states. Due to that non-unique absorbing phase the 
universality hypothesis of Janssen and Grassberger can not be applied. There- 
fore, the critical behavior of the pair contact process was intensely investigated 
by simulations (including [39,40,41,42,43]). It was shown numerically that the 
critical scaling behavior of the 1 + 1-dimensional pair contact process is charac- 
terized by the same critical exponents [31,40] as well as by the same universal 
scaling functions as directed percolation [26]. In particular the latter result pro- 
vides a convincing identification of the universal behavior. Thus despite the 
different structure of the absorbing phase, the 1 + 1-dimensional pair contact 
process belongs to the directed percolation universality class. This numerical 
evidence confirms a corresponding renormalization group analysis predicting 
DP universal behavior [44] in all dimensions. But the scaling behavior of the 
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PCP in higher dimension is still a matter of controversial discussions. A re- 
cently performed renormalization group analysis conjectures that the pair con- 
tact process exhibits a dynamical percolation-like scaling behavior [45,46]. A 
dynamical percolation cluster at criticality equals a fractal cluster of ordinary 
percolation on the same lattice. Thus, the dynamical percolation universal- 
ity class [47,48,49] differs from the directed percolation universality class. In 
particular the upper critical dimension equals D c = 6 instead of D c = 4 for 
DP. Furthermore, the dynamical scaling behavior of the PCP is a matter of 
controversial discussions. Deviations from the directed percolation behavior 
due to strong corrections to scaling [41] as well as non-universal behavior of 
the exponents [50] are observed in low dimension systems. 

So far, the investigations of the PCP are limited to the 1+1-dimensional [26,39,40,41,42,43] 
and 2 + 1-dimensional [51] systems. In this work we consider for the first time 
the PCP in higher dimensions and identify the scaling behavior via universal 
scaling functions. 

2.4 Threshold transfer process 

The threshold transfer process (TTP) was introduced in [32]. Here, lattice 
sites may be empty (n — 0), occupied by one particle (n — 1), or occupied by 
two particles (n = 2). Double occupied lattice sites are considered as active. 
In that case both particles may move to the left (1) or right (r) neighbor if 
possible, i.e., 

ri] — > rii + 1 if n-i < 2 , 

n v — > n T + 1 if n r < 2 . (9) 

Additionally to the particle movement, creation and annihilation processes are 
incorporated. A particle is created at an empty lattice site (0 — > 1) with prob- 
ability r whereas a particle annihilation (1 — > 0) takes place with probability 
1 — r. In the absence of double occupied sites the dynamics is characterized by 
a fluctuating steady state with a density r of single occupied sites. The density 
of double occupied sites is identified as the order parameter of the process, and 
any configuration devoid of double occupied sites is absorbing. The probabil- 
ity r controls the particle density, and a non-zero density of active sites occurs 
only for sufficiently large values of r. In contrast to the infinitely many frozen 
absorbing configurations of the pair contact process, the threshold transfer 
process is characterized by fluctuating absorbing states. Nevertheless steady 
state numerical simulations of the 1 + 1-dimensional threshold transfer pro- 
cess yield critical exponents that are in agreement with the corresponding DP 
values [32]. So far, no systematic analysis of the TTP in higher dimensions 
was performed. 
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In this work we limit our investigations to the 2 + 1-dimensional TTP. Anal- 
ogous to the 1 + 1-dimensional case, both particles of a given active site are 
tried to transfer to randomly chosen empty or single occupied nearest neigh- 
bor sites. Furthermore, we apply an external field that is conjugated to the 
order parameter. In contrast to the models discussed above the conjugated 
field can not be implemented by particle creation. Particle creation with rate 
h would affect the particle density, i.e., the control parameter of the phase 
transition. But the conjugated field has to be independent of the control pa- 
rameter. Therefore, we implement the conjugated field by a diffusion-like field 
that acts by particle movements. A particle on a given lattice site moves to 
a randomly selected neighbor with probability h, if n < 2. Thus the conju- 
gated field of the TTP differs from the conjugated field of the Domany-Kinzel 
automaton, the contact process, and the pair contact process. 

2.5 Ziff-Gulari- Bar shad model 

Another model exhibiting a directed percolation-like absorbing phase tran- 
sition is the Ziff-Gulari-Barshad (ZGB) model [33]. This model mimics the 
heterogeneous catalysis of the carbon monoxide oxidation 

2 CO + 2 — ► 2C0 2 (10) 

on a catalytic material, e.g. platinum. The catalytic surface is represented 
by a square lattice where CO or 2 can be adsorbed from a gas phase with 
concentration y for carbon monoxide and 1 — y for oxygen, respectively. The 
concentration y is the control parameter of the model determining the den- 
sity of the components on the catalytic surface. Adsorbed oxygen molecules 
dissociate at the catalytic surface into pairs of O atoms. It is assumed that 
the lattice sites are either empty, occupied by a CO molecule, or occupied 
by an O atom. Adjacent CO and O react instantaneously and the resulting 
C0 2 molecule leaves the system. Obviously, the system is trapped in absorb- 
ing configurations if the lattice is completely covered by carbon monoxide or 
completely covered by oxygen. The dynamics of the system is attracted by 
these absorbing configurations for sufficiently large CO concentrations and for 
sufficiently large 2 concentrations. Numerical simulations show that catalytic 
activity occurs in the range 0.390 < y < 0.525 [52] only. The system under- 
goes a second order phase transition to the oxygen passivated phase whereas 
a first order phase transition takes place if the CO passivated phase is ap- 
proached. In particular, the continuous phase transition is expected to belong 
to the universality class of directed percolation [53]. This conjecture is sup- 
ported by numerical determinations of certain critical exponents [52,54]. At 
first glance, it might be surprising that the ZGB model exhibits directed per- 
colation like behavior since the ZGB model is characterized by two distinct 
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chemical components, CO and O. But the catalytic activity is connected to 
the density of vacant sites, i.e., to a single component order parameter [53]. 
Thus the observed directed percolation exponents are in full agreement with 
the universality hypothesis of Janssen and Grassberger. But one has to stress 
that the ZGB model is an oversimplified representation of the catalytic carbon 
monoxide oxidation. A more realistic modeling has to incorporate for example 
mobility and desorption processes as well as a repulsive interaction between 
adsorbed oxygen molecules (see e.g. [33,55]). These features affect the critical 
behavior and drive the experimental system out of the directed percolation 
universality class. 

In this work, we focus on the two-dimensional ZGB model and determine cer- 
tain scaling functions for the first time. Therefore, we apply an external field 
conjugated to the order parameter. Since the order parameter is connected to 
the density of vacant sites of the catalytic reaction, the conjugated field could 
be implemented via a desorption rate h of adsorbed oxygen molecules. 

In summary, we investigate the scaling behavior of five different models span- 
ning a broad range of interaction details, such as different update schemes (ran- 
dom sequential as well as parallel update), different lattice structures (simple 
cubic and bcc lattice types), different inactive backgrounds (trivial, fluctuat- 
ing or quenched background of inactive particles), different structures of the 
absorbing phase (unique absorbing state or infinitely many absorbing states), 
as well as different implementations of the conjugated field (implemented via 
particle creation, particle diffusion or particle desorption). Nevertheless we 
will see that all models are characterized by the same scaling behavior, i.e., 
they belong to the same universality class. 



3 Steady state scaling behavior 

In this section we consider the steady state scaling behavior close to the 
transition point. Therefore, we performed steady state simulations of the five 
models described above. In particular, we consider the density of active sites 
p a = (L~ D N a ), i.e., the order parameter as a function of the control param- 
eter and of the conjugated field. Analogous to equilibrium phase transitions, 
the conjugated field results in a rounding of the zero- field curves and the or- 
der parameter behaves smoothly as a function of the control parameter for 
finite field values (see figure 1). For /i^Owe recover the non-analytical order 
parameter behavior. Additionally to the order parameter, we investigate the 
order parameter fluctuations Ap a = L D ((pl) — (p a ) 2 ) and its susceptibility x- 
The susceptibility is obtained by performing the numerical derivative of the 
order parameter p a with respect to the conjugated field (8). Similar to the 
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Fig. 1. The order parameter /? a and its fluctuations Ap a (inset) of the 
2 + 1-dimensional transfer threshold process (TTP) on a square lattice for various 
values of the field (from h = 10 -5 to h = 10~ 3 ). The symbols mark the zero-field 
behavior. The data are obtained from simulations on various system sizes with pe- 
riodic boundary conditions. 

equilibrium phase transitions, the fluctuations and the susceptibility display a 
characteristic peak at the critical point. In the limit h — > this peak diverges, 
signalling the critical point. 

In the following, we take only those simulation data into account where the 
(spatial) correlation length £j_ is small compared to the system size L. In 
that case, the order parameter, its fluctuations, as well as the order parame- 
ter susceptibility can be described by the following generalized homogeneous 
functions 



with the order parameter exponent /3, the field exponent a (corresponding 
to the gap exponent in equilibrium), the fluctuation exponent 7', and the 
susceptibility exponent 7. Here, h denotes the conjugated field and 5p denotes 
the distance to the critical point, e.g. 5p = (A — A c )/A c for the contact process, 
5p = (p — p c )/p c for site-directed percolation and for the pair contact process, 
5p = (r — r c )/r c for the threshold transfer process, etc.. The so-called non- 
universal metric factors a p , a h , a A and a x contain all non-universal system 
dependent features [56] (e.g. the lattice structure, the range of interaction, 
the used update scheme, as long as the interaction decreases sufficient rapidly 




(11) 
(12) 

(13) 



a A Ap a (5p, h) 



A 7 D(a p 5p \,a h h\°~) , 
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as a function of separation, etc.). Once the non- universal metric factors are 
chosen in a specified way (see below), the universal scaling functions R, D, 
and X are the same for all systems within a given universality class. The above 
scaling forms are valid for D ^ D c . At the upper critical dimension D c they 
have to be modified by logarithmic corrections [57,18]. 

Throughout this work we norm the universal scaling functions by R(l, 0) = 1, 
R(0, 1) = 1, and D(0, 1) = 1. In that way, the non-universal metric factors a p , 
a h , and a A are determined by the amplitudes of the power-laws 



p*(6p,h = 0) ~ (a p Spf, (14) 
p>(6p = 0,h) ~ {a h hf' a , (15) 
a A Ap^5p = 0, h) ~ (a, h)~^ >" . (16) 

Taking into consideration that the susceptibility is defined as the derivative 
of the order parameter with respect to the conjugated field (8) we find 

X(x,y) = d y R(x,y), a x = a' 1 , (17) 



as well as the scaling law 




Fig. 2. The universal scaling function R(x, 1) of the directed percolation universal- 
ity class. The data are plotted according to equation (21). All models considered are 
characterized by the same universal scaling function, an impressive manifestation 
of the robustness of the directed percolation universality class with respect to vari- 
ations of the microscopic interactions. Neglecting the non-universal metric factors 
a p and a h each model is characterized by its own scaling function (see inset). For all 
models the scaling plots contain at least four different curves corresponding to four 
different field values (see e.g. figure 1). The circles mark the condition R(0, 1) = 1. 
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a p 5p (a h h)' 1 ° 

Fig. 3. The universal scaling function 1) of the directed percolation universality 
class in various dimensions. The two- and three-dimensional data are vertically 
shifted by a factor 8 and 64 in order to avoid overlaps. The circles mark the condition 
£(0,1) = 1. 

This scaling law corresponds to the well known Widom law of equilibrium 
phase transitions. Furthermore, comparing equation (13) for 5p = to the 
definition of the susceptibility 

a xX (Sp,h) ~ (aJi)-^X(0,l), x = d h p> = d h (a h h)^ (19) 
leads to 

X(0,1) = A. (20) 

This result offers a useful consistency check of the numerical estimates of the 
susceptibility. Furthermore, it is worth mentioning that the validity of the 
scaling form (11) implies the required singularity of the susceptibility (8), i.e., 
it confirms that the applied external field is conjugated to the order parameter. 

Choosing a h h\ a = 1 in equations (11-13) we obtain the scaling forms 



p^5p, h) ~ (a h hf' a R(a p 5p (a h h)' 1 ^, 1) , (21) 

a A Ap^5p, h) ~ K h)-^'" D(a p 5p (a h h)' 1 ^ , 1) , (22) 

«x xfe /i) ~ (a, h)-~</° X(a p 5p (a h h)' 1 ^, 1) . (23) 
Thus plotting the rescaled quantities 

Pa (a h h)~^ , a A Ap c (a h hy''° , a x X (a h h)^'" (24) 
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Fig. 4. The universal scaling functions of the susceptibility X(x, 1) and of the 
fluctuations D(x, 1) (inset) in various dimensions. In case of the susceptibility the 
two- and three-dimensional data are vertically shifted by a factor 3/2 and 9/4 in 
order to avoid overlaps. The circles mark the condition X(0, 1) = (3 ju and reflect 
the accuracy of the performed analysis. The two- and three-dimensional fluctuation 
data are vertically shifted by a factor 8 and 64 in order to avoid overlaps. The circles 
mark the condition D(0, 1) = 1. 

as a function of the rescaled control parameter x = a p 5p (a h h)~ l l (T the data of 
all systems belonging to the same universality class have to collapse onto the 
universal curves R(x, 1), D(x, 1), and X(x, 1). This is shown for R(x, 1) in fig- 
ure 2 where the rescaled order parameter is plotted for various two-dimensional 
models. For this purpose the best known estimates for the critical exponents, 
as given in Table 1, are used. As can be seen the data of all considered mod- 
els collapse onto the same scaling function, clearly supporting the assumption 
that all models belong to the universality class of directed percolation. Fur- 
thermore, the data-collapse confirms the accuracy of the numerically estimated 
values [54] of the critical exponents. 

To limit the numerical effort we skip in the following analysis the TTP and 
ZGB and focus to the scaling behavior of the CP, sDP, and the PCP in various 
dimension. The universal scaling function R(x, 1) is displayed in figure 3 for 
D = 1,2,3. For each dimension, the data of the three models collapse onto 
the unique scaling function. As expected the scaling functions vary with the 
spatial dimensionality below for D < D c . 

The universal scaling functions of the order parameter fluctuations and the 
order parameter susceptibility are shown in figure 4. The susceptibility is ob- 
tained by performing the numerical derivative of the order parameter with 
respect to the conjugated field. The perfect data-collapses confirm the scaling 
forms (12,13). All scaling functions exhibit for D — 1,2,3 a clear maximum 
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signalling the divergence of Ap a and \ & t the critical point. The susceptibility 
data fulfill equation (20), reflecting the accuracy of the performed analysis. 
In summary, all considered models are characterized by the same critical ex- 
ponents and the same steady state scaling functions R, D, and X. Thus the 
steady state scaling behavior of the contact process, the Domany-Kinzel au- 
tomaton, the pair contact process, the threshold transfer process, as well as 
of the Ziff-Gulari-Barshad model belong to the same universality class for 
D = 1,2,3. 



3. 1 Mean field regime 



Mean field theories of all models considered can be simply derived and are 
well established (see for example [58,59,60,34,61,26]). These mean field the- 
ories present not only some insight into the critical behavior, they become 
valid above the upper critical dimension D c . Thus analytical expressions for 
the scaling functions become available for D > D c . In the case of directed 
percolation the mean-field scaling functions are given by (see e.g. [7,19]) 



x I /x x 2 



R**(x,y) = - + p + (^-J , (25) 

D MF (x,y)= , (26) 

\/y + (z/2) 2 

X MF (x,y) = 1 (27) 

2 sjy + (:r/2) 2 

i.e., the mean-field exponents are fi MF = 1, a MF = 2, 7 MF = 1, and 7^ F = 0. 

The scaling behavior of the fluctuations deserves comment. The exponent 

7^ IF = corresponds to a jump of the fluctuations and the scaling form (12) 
reduces to 

a A Ap,(5p, h) ~ D MF (a p 5p A, a h h A CT ) . (28) 



Using again D MF (0, 1) = 1, the non-universal metric factor a A is determined 

by 

°* = A A (*p=o,Ar (29) 



Numerical data of the five-dimensional models are presented in figure 5. A 
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Fig. 5. The universal scaling function of the order parameter R(x, 1) and the fluc- 
tuations D(x, 1) (inset) for D = 5. The five-dimensional data agree with the corre- 
sponding mean field functions (30) and (31), respectively. 

good data-collapse of the rescaled data [rescaled according to equations (21)- 
(23)] with the universal mean field scaling functions 




(30) 

D MF (x, 1) = 1 + (31) 

2 y/l + (x/2) 2 



is obtained. 



The observed agreement of the data of the five- dimensional models with the 
mean field scaling functions R MF (x, 1) and D MF (x, 1) leads to the result that 
the upper critical dimension is less than five, in agreement with the renor- 
malization group conjecture D c = 4 [62,63]. This is a non-trivial result. The 
value of D c is often predicted by field theoretical treatments, and a reliable 
and independent determination of the upper critical dimension is therefore of 
particular interest. For example, two contrary field theories conjecture for the 
pair contact process D c = 4 [44] and D c = 6 [45,46], respectively. The latter 
result is in clear contradiction to the observed mean field behavior in D = 5. 



3.2 Upper critical dimension 



The scaling behavior at the upper critical dimension is characterized by mean 
field power-laws modified by logarithmic corrections. Recent numerical inves- 
tigations [64,57,19,65] as well as analytical results [18] reveal that the concept 
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of universal scaling functions can also be applied to the upper critical dimen- 
sion. For example, the order parameter is expected to obey the scaling form 
(all terms in leading order) 

a aP ^5p,h) ~ A^ MF |lnA| A R{a p 5p A | In A| n , a h h A<™ F | In A| H ) . (32) 



Greek capitals will be used in the following to denote the logarithmic correction 
exponents. According to this scaling form, the order parameter at zero field 
(h = 0) and at the critical density (Sp = 0) is given by 



a a p a (<5p, h = 0) ~a p 5p \ \na p 5p\ B #(1,0), 



a,p,(5p = 0,h)~^aji In ^aji 



£(0,1), 



(33) 
(34) 



with B = II + A and E = H/2 + A. Similar to D ^ D c the normalization 
R(0, 1) = R(1,0) = 1 is used. Furthermore, the scaling behavior of the equa- 
tion of state is given in leading order by 



a a p a (<5p, h) ~ Ja h h 



In J a h h 



R(x, 1) , 



(35) 



where x denotes the scaling argument 



x = a p 5pJa h h 



In J a h h 



(36) 



with^ = II-H/2 = B-£. 

In case of directed percolation it is possible to confirm the scaling form (32) by 
a renormalization group analysis [18], yielding in addition the values A = 7/12, 
II = —1/4, and H = —1/2. Thus the scaling behavior of the equation of state 
is determined by B = E = 1/3 and ^ = 0. 

Similarly to the order parameter the following scaling form is used for the 
fluctuations [57,19] 

a A Ap a (5p, h) ~ A 7 mf | In A| K D(a p 5p A | In A| n , a h h \~ a \ In A| H ) . (37) 



Taking into account that numerical simulations show that fluctuations remain 
finite at the critical point (i.e. 7^ IF = and K = [19]) the scaling form 

a A Ap^5p,h) ~ D(x,l) (38) 
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Fig. 6. The universal scaling functions of the order parameter and its fluctuations 
(inset) at the upper critical dimension D c = 4. The logarithmic correction exponents 
are given byB = £ = l/3 [18] and ^ = 0. For both models considered, the scaling 
plots contain at least three different curves corresponding to three different field 
values. The circles mark the condition .R(0, 1) = 1 and D(0, 1) = 1, respectively. 



is obtained, where the scaling argument x is given by equation (36) with ^ = 
0. The non-universal metric factor a A is determined again by the condition 
£(0,1) = 1. 

In that way, the scaling behavior of the order parameter and its fluctuations 
at D c is characterized by two exponents (B = 1/3 and E = 1/3) and four 
unknown non-universal metric factors (a s ,,a p ,a h ,a A ). The detailed analysis 
for sDP data is described in [19]. Here we present a corresponding universal 
scaling plot containing data of two models, namely sDP and the CP. The 
corresponding scaling plots are displayed in Fig. 6 and show that the concept 
of universal scaling functions can be applied to the upper critical dimension. 

Notice that no data-collapse is obtained if logarithmic corrections are ne- 
glected, i.e., for B = E = 0. Thus, at least the leading logarithmic corrections 
have to be taken into account in order to study steady state scaling functions. 
It is therefore remarkable that recently performed off-lattice simulations of 
the dynamical scaling behavior at D c = 4 reveal that logarithmic correc- 
tions of higher orders (e.g. O(lnlnt)) are necessary to describe the numerical 
data [66]. Although the steady state results presented here are quite convinc- 
ing, we expect that even better results are obtained by incorporation higher 
order corrections. 
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4 Dynamical scaling behavior 



In this section we discuss the dynamical scaling behavior close to the transi- 
tion point. We limit our attention to the so-called activity spreading, generated 
from a single active seed. Starting with the seminal work of [67], measurements 
of activity spreading have been widely applied in the last two decades. In par- 
ticular, they provide very accurate estimates of the critical values p c , A c ,. . . as 
well as of the exponents 5, 9, and z (see below). Here, we will focus on the 
scaling functions of the survival probability P a and of the average number of 
active sites iV a . At criticality both quantities obey the power-laws 

a P P a ~ (a t t)~ 5 , a N N^ ~ (a t t) e , (39) 

where is termed the critical initial slip exponent [68]. Sufficiently close to 
the critical point, P a and iV a are expected to obey the scaling forms 

a P P a (<5p, L, t) ~ A"'"" P pb » A, a L L , a t t A""" ) , (40) 
a N N a (5p, L, t) ~ \ ev w N phc (a P 5p A, a L L , a t t A""" ) , (41) 

where the time t and the system size L are incorporated as additional scaling 
fields. The scaling power of t has to equal the scaling power of the correlation 
time £|| oc 8p~ v w, whereas the scaling power of the system size is given by 
the spatial correlation length exponent u± (£j_ oc 5p~ u± ). In contrast to the 
previous section, the system size L has to be taken into account because the 
power-law behaviors (39) are limited by the finite system size. The index 
pbc indicates that the universal finite-size scaling functions depend on the 
particular choice of the boundary conditions as well as on the system shape (see 
e.g. [69,70,71,72,73]). But different lattice structures are contained in the non- 
universal metric factors. 

Choosing a t t\~^ = 1 the power-laws (39) are recovered for P pbc (0, oo, 1) = 1 
as well as iV pbc (0, oo, 1) = 1. The finite-size scaling forms are obtained by 
setting a L L \~ u± = 1, yielding 

a P P a (0, L, t) ~ (a L L)~ Sz P pbc (0, 1, a t t (a L L)~ z ) , (42) 
a N iV a (0, L, t) ~ (a L L) ez A> pbc (0, 1, a t t (a L L)- z ) . (43) 

The scaling functions P pbc (0, 1, x) and A r pbc (0, 1, x) are expected to decay expo- 
nentially for t ^> t PSS whereas they exhibit an algebraic behavior for t t PSS , 
with t FSS = a t " 1 (a L L) Jy± . 
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Fig. 7. The universal scaling functions P pbc and iV pbc of activity spreading for 
various dimensions. In case of the pair contact process the simulations are started 
from a natural configuration of inactive particles. System sizes L = 64, 128, 256, 512 
are considered for D = 1, L = 64,128,256,512 for D = 2, and L = 16,32,64,128 
for D = 3. The dashed lines corresponds to the power-law behavior of the infinite 
system x~ s and x 9 , respectively. 

Performing activity spreading simulations of site directed percolation and of 
the contact process the initial seed is implemented by a single particle on 
an empty lattice. For absorbing phase transitions with non-trivial absorbing 
states, like the pair contact process, the scaling behavior depends upon the 
nature of the initial configuration [40]. In that case, spreading activity simu- 
lations have to be performed at the so-called natural density of inactive par- 
ticles [40,32,35]. Starting with a random configuration, an absorbing state at 
criticality is prepared by the dynamics. An active seed is created for example 
by adding or moving a particle. The resulting activity relaxation is monitored 
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Table 1 

The critical exponents of directed percolation for various dimensions D. In D = 1, 
the exponents 7, i/±, and i/y are obtained from a series expansion by Jensen [85]. 
For D = 2 and D = 3 activity spreading simulations are performed yielding S, 9, 
as well as z [54,86]. Additionally, the exponent i/m is determined [37,86] in order to 
estimate the full set of exponents via scaling laws. 





D = 1 [85] 


D = 2 [54,37] 


D = 3 [86] 


Mean field 





0.276486(8) 


0.5834 ± 0.0030 


0.813 ± 0.009 


1 




1.096854(4) 


0.7333 ± 0.0075 


0.584 ±0.005 


1/2 




1.733847(6) 


1.2950 ± 0.0060 


1.110 ±0.010 


1 


a 


2.554216(13) 


2.1782 ±0.0171 


2.049 ± 0.026 


2 


i 


0.543882(16) 


0.2998 ± 0.0162 


0.126 ±0.023 





7 


2.277730(5) 


1.5948 ± 0.0184 


1.237 ±0.023 


1 


5 


0.159464(6) 


0.4505 ± 0.0010 


0.732 ± 0.004 


1 


e 


0.313686(8) 


0.2295 ± 0.0010 


0.114 ±0.004 





z 


1.580745(10) 


1.7660 ± 0.0016 


1.901 ±0.005 


2 



until it ceases. Then we start a new measurement from a random configura- 
tion and so on. Therefore, the numerical effort is significantly increased for 
systems exhibiting non-trivial absorbing states and only small system sizes 
are available by simulations. But nevertheless convincing data-collapses, in- 
cluding the pair contact process data, are obtained and the corresponding 
universal scaling functions are shown in figure 7. The values of the exponents 
used are listed in table 1. In summary, activity spreading from a localized seed 
is characterized by the same universal scaling functions P pbc and N pbc for all 
considered models. Thus the CP, sDP, and especially the PCP display the 
same dynamical scaling behavior at criticality. 



5 Conclusions 

Similar to equilibrium critical phenomena, the great variety of non-equilibrium 
phase transitions can be grouped into different universality classes. Each non- 
equilibrium universality class is characterized by a certain symmetry which 
is often masked within the Langevin equation approach, but it is expressed 
clearly within the corresponding path integral formulation (see e.g. [2,15] and 
references therein). For example, non-equilibrium critical systems belong to 
the directed percolation universality class if the associated absorbing phase 
transition is described by a single component order parameter and if the cor- 
responding coarse grained system obeys the rapidity reversal symmetry (at 
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least asymptotically). This is clearly demonstrated in figure 2 where the data 
of five different models is shown. Although the models exhibit different in- 
teraction details their rescaled data collapse onto a unique universal scaling 
curve. This scaling plot is an impressive manifestation of the robustness of 
the DP universality class. On the other hand, it allows the identification of 
the irrelevant parameters, i.e., those parameters which do not affect the uni- 
versality class. As expected, the used update scheme, the lattice structure 
and the different implementation schemes of the conjugated field do not af- 
fect the scaling behavior. But it is worth mentioning that the structure of the 
absorbing phase was considered so far as a relevant parameter. For example, 
it was expected that models with infinitely many absorbing states belong to 
a different universality class than DP. The observed (steady state and dy- 
namical) DP-scaling behavior of the pair contact process at criticality reveals 
that the number of absorbing states does not contribute to the asymptotic 
scaling behavior at criticality. For the sake of completeness, we mention that 
different universality classes than DP occur if the rapidity reversal is broken, 
e.g. by quenched disorder [1,74,75,76,77,78,79], or additional symmetries [3,4] 
such as particle-hole symmetry (compact directed percolation [30,80]), or par- 
ticle conservation (Manna universality class [81,82,57]) or parity conservation 
(for example branching annihilating random walks with an even number of 
offsprings [83,84]). 

We would like to thank P. Grassberger for communicating his results prior to 
publication. Furthermore, we thank H.-K. Janssen and R. M. Ziff for instruc- 
tive discussions and useful comments on the manuscript. 
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